Melting of persistent double-stranded polymers 
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Motivated by recent DNA-pulling experiments, we revisit the Poland-Scheraga model of melting 
a double-stranded polymer. We include distinct bending rigidities for both the double-stranded 
segments, and the single-stranded segments forming a bubble. There is also bending stiffness at the 
branch points between the two segment types. The transfer matrix technique for single persistent 
chains is generalized to describe the branching bubbles. Properties of spherical harmonics are then 
exploited in truncating and numerically solving the resulting transfer matrix. This allows efficient 
computation of phase diagrams and force-extension curves (isotherms). While the main focus is on 
exposition of the transfer matrix technique, we provide general arguments for a reentrant melting 
transition in stiff double strands. Our theoretical approach can also be extended to study polymers 
with bubbles of any number of strands, with potential applications to molecules such as collagen. 
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I. INTRODUCTION 

Single-molecule micromanipulation techniques have 
opened up new opportunities for measurements and stud- 
ies of polymers. Smith et al. pioneered [l| stretching ex- 
periments of double-stranded DNA (dsDNA) and, along 
with others, observed that at high forces of about 65pN, 
DNA extends to 1.7 times its contour length 0, 0,11, 
These investigators believe that the stretching transforms 
B-DNA, which is DNA in its natural state, to a new, 
extended state, named S-DNA. Modeling studies and 
simulations were carried out to characterize this puta- 
tive new state of DNA. 0, H, 0] Subsequently, Storm and 
Nelson [l(| proposed a statistical model of DNA as a 
discrete persistent chain (DPC) with two monomer fla- 
vors of different lengths and stiffnesses, and fit their pa- 
rameters successfully to experimental data. However, 
Williams, Rouzina, Bloomfield, and co-workers have ar- 
gued on the basis of their own experiments that S-DNA 
is not a new state of the molecule, but merely DNA 
that is melted to two sing le-stranded DNA (ssDNA) 
fragments . [Til [l2l. fl3l. [l 5l. 1 1 6l| Furthermore, they deem 
the aforementioned modeling and simulations of S-DNA 
as contradicting experimental data. Furthering this con- 
troversy, Cocco et al.fl7j reexamine the experimental 
data and argue in favor of S-DNA, Whitelam et al.[l8| do 
so based on kinetics, while Piana [19| observes melting in 
simulations of short stretches of DNA. 

In 1966 Poland and Scheraga[2p|| introduced a simple 
statistical model for the melting of the dsDNA to two 
ssDNA fragments, which has proved quite illuminating. 
In this model, configurations of partially melted DNA 
are represented by alternating segments of dsDNA, and 
denatured pairs of single strands forming 'bubbles.' To 
make the model analytically tractable, certain features 
of DNA such as excluded volume, bending rigidity, and 



sequence inhomogeneity are typically left out. With the 
later inclusion of excluded volume effects, the model is 
well suited for characterizing the nature of the melting 
transition, and its universality. For comprehensive (but 



older) reviews see Rcfs. 
described in, e.g. Ref. 
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[22j; some newer results are 
. More recently, the phase 
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diagram of the model has been studied in the presence 
of a stretching force [H, [25[. This is important, since 
even the experiments disputing the formation of S-DNA 
at 65pN do observe melting induced stretching at other 
forces[5, 6|. The effect of bending rigidity is still left out 
in the newer studies, making comparisons to experiment 
questionable. The aim of this paper is to facilitate the 
ongoing debate by providing a model that accounts for 
the bending rigidity of the polymer (while leaving out 
excluded volume effects). 



While we hope that our results and phase diagrams 
provide an additional perspective into this system, our 
main accomplishment is the extension of the transfer ma- 
trix method used for a single persistent polymer (worm- 
like chain) to the melting of a double-stranded poly- 
mer. The remainder of the paper is an exposition of 
our method, and is organized as follows. The general- 
ized Poland-Schcraga model with three types of bending 
rigidity is introduced in Sec. Ill A[ and the correspond- 
ing three contributions to transfer matrices are developed 
in Sec. Ill Bl As described in Sec. IIII1 numerical results 
can be obtained by truncating the resulting transfer ma- 
trices in a basis of spherical harmonics. In particular, 
we provide phase diagrams (in force and temperature) 
and force-extension curves, along with the native (dou- 
ble stranded) fraction. We augment numerical results 
with physical explanations of the observed trends. In 
particular, we provide a rather general characterization 
of the slop of the phase boundary which explains the 
potential reentrant character of force induced melting. 
Various technical details of the calculation are relegated 
to the Appendices. 
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II. MODEL 
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FIG. 1: A typical polymer configuration of our model, as 
depicted here, consists of segments R, solid arrows, which we 
imagine to be dsDNA, alternating with 'bubbles' B made of 
two strands of ssDNA (light arrows). The two segment types 
have unit (monomer) lengths \f\ or \b\, and bending costs 
of Jr or Jb, respectively. There is an additional bending 
constant Jj, at the branching points, and a weight wj = e ej 
for each joint. The energetic advantage (binding energy) of 
the R segments is represented by a weight wr = e tR per step. 



As illustrated in Fig. [TJ a typical configuration of our 
model polymer consists of an alternating sequence of na- 
tive segments R, and locally molten pairs of strands form- 
ing a bubble B. Successive segments are indexed by i, 
and contain Nrj or 2Ns.i monom ers, respectively. In 
the original Poland-Scheraga model [20] , the R segments 
were treated as stiff 'r'ods. We treat these segments as 
scmi-flcxiblc chains, such that the energy of a segment of 
Nrj monomers is given by 



N R -1 

-f3E R = ^2 [ J R • fj+i +?j ■ f 

3=1 

+ rN R ■ f + N R e R . 



(1) 



Here, fj = \r\fj is the displacement of the j'th 'monomer' 
of the segment, all of which have equal length, but may 
point in any direction. The coupling J R parameterizes 
the cost of bending neighboring monomers. The force 
/ stretches the polymer, and e R is an additional contri- 
bution to the energy difference between a native R unit 
compared to the molted strands of B units. Note that 
(for each configuration, and discounting bending costs) 
the net energy difference between bound and unbound 
segments (the binding energy) is ksT(J R + e R ) per base- 
pair. (For ease of notation, the index i denoting the i'th 
R segment has been dropped from all variables above.) 

Similarly, the energy of a molten B region, described 



/ 



(2) 



N B -1 I 

-j3E' B = [JBb r b j+1 + b 

3 = 1 V 

N B -1 / 7 
3=1 \ / 

+ t>N B ■ 2 + %n b ■ 2 ■ 

Again, the implicit index i numbering the i'th B seg- 
ment has been omitted. The allowed configurations are 

constrained by Rb — Y^j=i bj = Ylj=i ^'ji to ensurc that 
the two branches of the bubble end at the same point. It 
is indeed this constraint (emphasized by the primed Eb) 
that allows distributing the energy cost of stretching by 
the force / symmetrically between the two branches. 

Finally, there is a joint when the iV/j^'th (last) element 
of the i'th R segment branches into the first elements of 
the i'th B segment, to which we associate an energy 



f3Ej tRB = Jj r Na ■ bi + J j f Nn ■ b[ + ej 



(3) 



Similarly at the point where the z'th B segment meets 
the (i + l)'th R segment, the energy is 



/3Ej,br = Jj b NB -fi + Jj V Nb ■ f i + ej 



(4) 



The overall energy of M alternating R-B segments of sizes 
{N R ,i,N B ,i} is thus 



(3E' [N RA ,N B ,i,N R , 2 ,--- ,7V, 



B,M\ 



M 

E 

i=l 



0E Rii + pEj >RB!i + f3E' B - + /3Ej. 



(5) 



BR.i 



(The above formula applies to configurations which start 
with an R segment and end with a B segment. We expect 
the results for long polymers to be independent of the 
choice of boundary conditions.) 

Computations are most easily performed in a grand 
canonical ensemble in which we sum over all possible 
polymer lengths, with a chemical potential fx//3 per 
monomer. The grand partition function is then calcu- 
lated from 



s 2 



N 



E ■ 



-PE'IN Ri1 ,-,N b ,m] 



(6) 

where N = Yli=i N R ,i + N B ,i is the native polymer 
length. The integrations are over all directions of the 
monomer vectors f, 6, and b' , provided that the bubble- 
closing constraints are satisfied. This can be ensured by 
inserting (^-functions for each bubble segment, as 



'N Bli JVb,« 



d 3 k 

(2^)3 



(7) 
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Transfer Matrix Formulation 



breviation a = (a, a') = ((l a , m a ), (l a i , m a i)), and 



The one-dimensional character of the energy in Eq. ([5|) 
suggests a transfer matrix approach to the problem. This 
is indeed a standard tool for the study of semi-flexible 
chains [l(| HH, H3, HI, H^] , but requires additional elabo- 
ration to treat the bubbles. Below, we shall develop step 
by step the contributions from the two segment types, 
and the joints in between, to the overall transfer matrix. 



^ e --- + J B b n - 1 -b n + J B b n -b n + 1 +b n -t+ib n -k+- 



S 2 (2tt) 3 d 2 b' n e --- + JBb' n _ 1 -b' n + J B ti n l' n + 1 +b' n -{-ib' n -k+- 



d 3 k 
(2^)3 



xYa(6 Il+ i) 



(10) 



where 



1. R segments 

The Boltzmann weight in Eq. involves a product 
of exponentials, similar in form to plane waves. Such 
exponentials can be expanded in a basis of spherical har- 
monics and Bessel functions, which then allows the inte- 
grations over the orientations f , b, and b' . For example, 
integrating over the unit vector f n of an R segment yields 



■ ■ + jR?n-l-'Tn + JRrn'rn+l+fn- f +CR+H+- ■ ■ ^2 » 

S 3 " (8) 

= [---W»-i)] i T R) a ,p [Yp(f n+1 )---} 7 

where summation over repeated indices is implied. Greek 
letters stand for elements of the angular momentum basis 
\l, m), e.g., a stands for (l a , m a ), and the transfer matrix 
elements are 

(T R ) a , = (47r) 2 C aJ ,XJ fi )i 7 (|/]H)F 7 *(/) w R z . (9) 

Here, i a is the modified spherical Bessel function of the 
first kind of order l a ; C a p a = f sa Y a (r)Yg(r)Yy(r) d 2 f 
is closely related to tabulated Gaunt coefficients, which 
can be expressed in terms of Wigner 3j-symbols, see Ap- 
pendix rAJ and each unit of an R segment carrier a fu- 
gacity z = e M , and the binding weight w B = e CR de- 
fined earlier. To make the notation uniform and simple, 
a bar placed over an index of C, e.g. C Q ^ „, indicates 
that the corresponding spherical harmonic under the in- 
tegral shall be complex conjugated. The repeated 7 in- 
dex implies a (finite) sum. The expression simplifies if 
the force / is chosen to point along the z direction, in 
which case (Tff) Q:( g oc 6 matTng . Note that the transfer 
matrix is asymmetric, as we have included ^(Jr), but 
not i a (Jii) to avoid double-counting. 



2. B segments 

A similar computation for the two bubble strands 
yields two transfer matrices. These must be combined 
into one matrix to be usable in the later steps. Thus, 
the basis elements \l,m), \l',m') are combined into one 
product basis element \l, m) <g> \l', m!) with one-letter ab- 



(T B (k)) &J = z / Y a (by-i+ fbrk YZ(b)(4Tr) l0 (J B )d 2 b 



x / Y a ,(b')e 



l\„b -t—ib ' -k 



Y; i {b'){ATT)i( 3l {J B )d 2 b' 



(11) 



This transfer matrix is, in general, very big. If the 
spherical harmonic basis elements are cut off at some 
! max for numerical evaluation, there are (Z max + l) 2 basis 
elements to consider since m £ {— I, ■ ■ ■ ,1} for each < 
I < imax- This means that there are ((Z m ax + 1) 2 ) 2 product 
basis elements, which is the number of rows or columns 
of the transfer matrix T B (k)\ With a little trick this big 
matrix can be reduced to block diagonal form with the 
biggest submatrix having size (Z max + l) 2 : 

Consider the first integral in Eq. (fTTj) and let v = ik+f/2 
indicate the vector in the exponent. One can rotate the 
complex vector v into the z-direction, such that 

(4tt) J Y a (b)e j; -%*(b)ip(J B )d 2 b = 

(4tt)P QiAI (^j Y»{b)e t ^Y:{b)i v {J B )d 2 i\ = (12) 

{^) 2 V C ^ (c^,MJB)h(W\\b\)Y;(z)) v~l , 

where £>(</>, 9, 0) is the rotation matrix (Wigner-D ma- 
trix or Wigner-D function) for quantum mechanical an- 
gular momentum states. The spherical coordinate angles 
9 (from the z-axis) and (f> (from the x-axis) are the an- 
gles by which the z-direction rotates into the ^-direction. 
The rotation matrices are usually parametrized by Eu- 
ler angles. The first Eulcr angle corresponds to <j>, the 
second to 9, and the third to zero. Since v is com- 
plex, the angles of rotation are complex and \v\, which 
is the usual Euclidean norm of v, is also complex, as 
|v| 2 = v ■ v = / 2 /4 - k 2 + ik ■ f. Note that ip(J B ) 
in the first line can be replaced by i u {J B ) because ip 
only depends on lp while T> only mixes states with the 



same total angular momenta, i.e. T)v,$ oc (5/,^, 



Al- 



though the dependence on k is now both in T> as well as 
in i 7 , the problem is computationally much simpler since 

C/J,P,7^-y (^) OC ^m^jTTljy- 
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3. Joints 

Similar manipulations lead to transfer matrices at the 
branching points of 



{Tj. R B) aj3 = {^fc a ^c BIjl ^{\fM)y;{f) 

x ip(Jj)i/3>(Jj) Wr wj z , 



(13) 



and 



(Tj,BR)a,p = (4^) 2 C a , a ,^C^ 7 ^ 7 (|/||r|)y 7 *(/) 

(Jj)ia'(Jj) 



1 a (JB)ia'(JB) 



ip(jR) Wr Wj Z . 



(14) 



Once more the lack of symmetry between the two cases 
reflects our choice of including the bending energy from 
the next (but not the previous) segment in the transfer 
matrix. The absence of the joint energy Jj greatly sim- 
plifies the problem, as the R and B segments can then be 
treated independently. This happens because i a (0) = 
unless l a = 0. 

Combining these expressions, we can express the parti- 
tion function T in Eq. (J5]) in terms of the transfer matrices 
(the irrelevant prcfactors have been omitted), as 



M=l 



i=0 

^ rTj.rbFbTj.br 





M 


j Tj.BR. 





i - t r t 



R 1 J,RB^ bTjrR 



0.0 



where 



and 



r« = (i - t r ) 



d 3 k T B (k) 2 



(2tt) 3 l-T B (k) 



0.0 



(15) 



(16) 



(17) 



III. RESULTS 



In the grand canonical ensemble the average length of 
the polymer is given by (N) = <9 M log(r) = z<9 z log(r). 
We are interested in the limit of very long polymers, 
where (N) — > oo. This limit is obtained for a specific 
choice of the fugacity z = e p , such that: 

1. There are infinitely many repetitions of native 
(R) and molten (B) segments. This occurs for 
a value of z such that the largest eigenvalue of 
FrT,j.rb^bTj.br]{z) equals 1. 

2. The size of an individual bubble diverges. In this 
case the singularity arises from d z T B (z). 



3. The hypothetical possibility of an infinitely long 
native (R) segment does not arise, as 8 z Tr(z) only 
diverges for values of z that already cause one of 
the previous two cases to occur. 

As in the Poland-Scheraga model [2(|, case (1) corre- 
sponds to a partially melted double strand (mixed phase 
of R and B segments), while case (2) corresponds to 
a fully melted state comprised of one bubble (bubble 
phase) . 



A. Phase Diagram 

A phase transition between the two phases occurs when 
there are both infinitely many repetitions of R and B seg- 
ments, and the average bubble size diverges. Since our 
model has 6 parameters, we have to select appropriate 
subspaces for the display of phase diagrams. We choose 
to regard the bending energies Jr, Jr, Jj, and the joint 
Boltzmann weight wj — e ej as parameters, and display 
the phase diagrams as a function of the force / and the 
dimensionless energy £r = \og{wR). The latter may be 
regarded as a stand-in for an inverse temperature, since 
it is related to an actual energy after division by fc^T. 
(As discussed following Eq. ([I}, negative values of £_r(/) 
do not pose a problem, as the actual binding energy is 
k B T(J R + e R ).) 

Some typical phase boundaries c_r(/) are presented in 
Fig. [21 We used the Mathematica software package on 
an Intel Pentium 3 GHz desktop computer to obtain the 
phase diagrams, each of which took a few hours of com- 
putational time when we included partial waves up to 
I = 1 in the bubble (B) partition functions, and I = 5 in 
the native (R) partition functions. It should be possible 
to reduce the computational time significantly by using 
more appropriate software. The truncation of the trans- 
fer matrices at these partial waves was justified consider- 
ing the small bending parameters chosen for the bubbles 
and the joints. 

Each solid curve depicts the phase boundary for a par- 
ticular choice of parameters. All curves correspond to 
rather stiff R-scgmcnts with Jr = 5.0, but for different 
choices of the bending bending parameters Jr and Jj. 
The upper portion of each figure corresponds to the par- 
tially melted native state, which contains both R and 
B segments. There is no native R segment left in the 
lower portion, and the polymer is a single bubble below 
the phase transition line. The lowest curve in the top fig- 
ure (red in the online version) corresponds to no bending 
in the bubble or vertex; J B = J j = 0. As we increase Jr 
to larger values of 0.05 and 0.1 in the top figure (green 
and cyan, online), the bubble phase becomes more sta- 
ble. (The phase boundary moves up as indicated by the 
arrow.) If we now fix J B = 0.05, and increase J j to val- 
ues of 0.05 and 0.1 in the bottom figure (blue and brown 
online) we find that moves down, as indicated by 

the arrow. The mixed phase is stabilized by stiffening 
the joints between R- and B-segmcnts. 
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FIG. 2: The phase boundaries £«(/) separate the partially 
melted R-B phase (above the line) from the fully melted bub- 
ble phase (below the line) , for various combinations of Jb and 
J. j. Top: The bubble stiffness Jb is increased. Bottom: The 
joint stiffness Jj is increased. In both plots, the arrow indi- 
cates the direction in which the phase boundary moves. In 
all curves \r\ = 1.0, \b\ = 1.7, J R = 5.0, and ej = 1.0. 



Several features of these phase diagrams are now com- 
mented on in more detail. 



1. Explanation of the trends in phase diagrams 

With our choice of parameters, increasing the value 
of any stiffness parameter J, makes the corresponding 
segment (or joint) more favorable. This is because the 
corresponding Boltzmann weights are monotonically in- 
creasing functions of J (as are the modified spherical 
Bessel functions i a (J)). For example, for a bubble seg- 
ment the overall weight increases with Jb, despite the 
fact that there are fewer configurations (and hence re- 
duced) entropy for the stiffened and stretched bubble. 



These trends are further magnified at larger force as the 
stretched segments gain even more weight by aligning to 
the force, as can be easily seen in Fig. [5] 



2. Reentrance & the phase boundary at low and high forces 

An interesting feature of the phase boundary in Fig. [5] 
is its reentrance, namely for certain choices of en (—0.8 < 
en < —0.2) the mixed R-B phase is stable at intermediate 
values of the force, but melts at both weak and strong 
force. This reentrance is also present in another model of 
denaturation which incorporates excluded volume effects 
in the bubbles, but no bending rigidity [2~i| . 

This feature can be explained by examining the limit- 
ing behaviors of the phase boundary at small and large /. 
For large /, the polymer (whether in native or denatured 
state) is stretched along the direction of the force. The 
contribution from entropy is relatively small in this limit, 
and one can estimate the location of the phase boundary 
by comparison of energies: The energy of a fully stretched 
rod segment is Jr + e R + f\r\ per base-pair. If the two 
strands are separated the energy changes to 2Jb + f\b\. 
The transition occurs for en as 2Jb — Jr + (\b\ — |^)/, 
which has a positive slope since B strands have longer 
monomers and are favored by the force. 

As shown in Appendix [Bl this argument can be made 
more rigorous and extended to all cases where the con- 
tribution from the joints can be ignored. In such cases 
the slope of the phase boundary is given by dfe R . = 
(Lb)/ (Ng) - (L R )/(N R ), where denotes the av- 

erage length per monomer, calculated for each segment 
type (B segment or R segment) treated separately. 

At zero force, the average end-to-end extension of each 
segment, (Lx), is zero by symmetry. The extension for 
small / is linear, with a force constant (susceptibility) 
that is easily related to the variance of the end-to-end 
extension at / = 0. Since the change in free energy is pro- 
portional to / 2 , the phase boundary is also quadratic in 
this limit. In the absence of joint stiffness, the curvature 
of the transition line can be related to the difference in 
susceptibilities by d 2 f e R = (L%) C /(N B )-(L R ) C /(N R ), see 
Appendix [Bl for a derivation. Here, (L x ) c /(Nx) denotes 
the variance in length of rods or bubbles per monomer, 
computed for one rod segment or one bubble segment 
subject to the same fugacity and force as for the whole 
molecule. Since it is easier to rotate and align the more 
rigid R segments in the direction of the force, their gain 
from the force is larger, and small force favors the native 
double-stranded phase. 

The reentrance in the phase diagram of reference [24| , 
mentioned in the first paragraph of this subsection, can 
be explained with these expressions as well. In this pa- 
per different behaviors are obtained as a function of a 
parameter A, which determines how statistically favored 
joints are. They observe a reentrant phase diagram for 
A = 0.01 (disfavoring joints) but not for A = 1 (many 
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joints). The variance per monomer in the lengths of the 
bubbles is roughly constant ({L|} c oc (N B ) as in a ran- 
dom walk), but for stiff rods the variance grows as the 
average size {(L 2 R ) C oc (Nr) 2 as in a directed walk). For 
small A, there are few joints and longer rods just after the 
phase transition into the mixed bubble-rod phase, mak- 
ing the variance per monomer large. According to the 
above formula for B 2 £r this leads to reentrance. Note 
that excluded volume effects do not substantially modify 
this argument. 

Williams et al. experimentally observe (Figure 5 of 
Ref. a reentrance in the phase boundary when they 
fit their data to a simple model. Unfortunately, the area 
of interest is merely extrapolated and the transition is 
not probed at high enough temperatures and low forces 
(« 90°) to unambiguously verify reentrance. 



B. Force— Extension Isotherms 

An important probe of phase behavior comes from the 
force-extension curves, in which the end-to-end distance 
of the polymer is measured as a function of increasing 
force. Without loss of generality and to simplify cal- 
culations, these curves are obtained for Jb = Jj = 
(with Jji ~ 5.0, ej = 1.0). For comparison, we also 
plot the curves corresponding to the pure worm-like chain 
(WLC) model in black (containing only R segments for 
wj = 0). The plotted 'extension' is the average length 
along the force direction, made intensive and dimension- 
less by dividing by the number of monomers (N), and 
the monomer length |r| of the R segment, i.e. 



X = 



(N\r\) 



d f T 



(18) 



The two panels in Fig. [3] were selected to correspond 
to parameters with (top) and without (bottom) reen- 
trant melting. (Consider horizontal lines in Fig. [2] for 
€r = —0.4 and €r = 1.0, respectively.) In both cases, 
the force-extension curves for the double-stranded poly- 
mer track the behavior of the worm-like chain closely in 
the mixed R-B phase, but deviate significantly in the 
denatured phases; most pronouncedly for the re-entrant 
transition. 

As mentioned earlier, current experiments indicate 
that the WLC model describes the extension of DNA ac- 
curately for small applied forces, but fails at large forces 
due to the appearance of an over-stretched region. In 
Fig. [4] we probe the corresponding region in more detail 
for our model, exploring the effect of bending rigidities 
(for er = 1.0 with a single transition). The top panel 
depicts the effect of increasing the bubble stiffness Jb, 
which makes the transition region appear sharper. In- 
creasing the joint stiffness Jj (bottom panel) has the op- 
posite effect. An interesting feature of the bottom panel 
is that the trends in X(f) are not monotonic in Jj, de- 
creasing the extension for weaker force, and increasing it 
for larger force, leading to a crossing point in between. 
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FIG. 3: Comparison of the extension curves X(f) = 
(L z ) / (N\r\) (solid lines) for a double-stranded polymer, and 
a single stranded wormlike chain. The top panel corre- 
sponds to a case where melting is reentrant (e_R = —0.4), 
while there is a single denaturation transition in the bot- 
tom panel for en — 1.0 (cf. Fig. [2J. The fraction of na- 
tive (double stranded) polymer, O = {N r )/{N), is indicated 
by dashed lines. All curves correspond to Jb = Jj = 0, 
\r\ = 1.0, \b\ = 1.7, Jr = 5.0, and ej = 1.0 for the double 
strands. 



The reader should note that we have taken ej = 1.0 in all 
curves, making the joints favorable and common. This 
choice is made to exaggerate the effect of the joint bend- 
ing for display, as well as to broaden the phase transition 
in Fig. 01 thus highlighting the features of our model. A 
more realistic value, namely ej small or negative, gives 
qualitatively similar results. 



C. 9: Native (R) fraction 

Figure [3] also includes the native fraction 8 as a func- 
tion of /, depicted by the dashed curves. This is defined 
as the fractional amount of R segments in the polymer, 
which can be computed from 



9 



(Nr) 
(N) 



zd z T 



(19) 



Note that goes to zero continuously on approaching 
the bubble phase, underscoring the second order nature 
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FIG. 4: Detailed view of the dimensionless extension X — 
(L z ) I (N\f\) as a function of the force /, close to the phase 
transition, for various combinations of Jb and Jj. Top - 
The bubble stiffness Jb is increased. Bottom - The joint 
stiffness Jj is increased. In both plots, the arrows indicate 
the direction the extension curve moves. In all curves |r| = 
1.0, |6| = 1.7, J R = 5.0, and ej = 1.0. 



of the phase transition. It approaches zero rapidly, but 
in a linear fashion. This is because the phase transition 
in our model belongs to the same universality class as 
the classic Poland-Scheraga model [2(J. The addition of 
bending rigidity is irrelevant close to the phase transi- 
tion, and excluded volume effects (which do modify the 
universality [23l. |24|) arc not included in our model. 



the latter, as longer persistent segments are less suscep- 
tible to fluctuations and excluded-volume effects. It is 
thus necessary that comparison of models to experiment 
should include the effect of rigidity, as we have attempted 
here. More generally, our formalism can be extended to 
decribe the unravelling of any number of strands, for ex- 
ample from 1 to 3 in the case of collagen [3l|, [HJ . 
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APPENDIX A: GAUNT COEFFICIENTS 

The Gaunt coefficients are defined as 

C a ,0, y = [ Y a {f)Y p (f)Y 7 {f)d 2 f , (Al) 
Js 2 

where Y a is the spherical harmonic with indices (l a , m a ). 
If a bar is put on top of an index of C, the corresponding 
spherical harmonic in the integrand is complex conju- 
gated. The relation 

Y* m = (-irn_ m , (A2) 

can be used to relate a modified Gaunt coefficient with 
barred indices to one without barred indices. A well- 
known expression for the Gaunt coefficient in terms of 
Wigner 3j-symbols is [33| 

/ (2/ a + I)(2^ + l)(2/ 7 + I) 
" V 4^ (A3) 

f la l[3 ly\ I la 1/3 l-y \ 

y0 OJ \m a mp m 1 J 

Using the properties of the Wigner 3j-symbols one can 
restrict and simplify the sums appearing in the partition 
function. 



IV. DISCUSSION 

We have introduced a formalism to address the role of 
bending rigidity in the denaturation of double-stranded 
polymers, DNA providing a prime example. There has 
been some controversy on interpreting experimental re- 
sults for melting of DNA, or its denaturation by force. 
There is strong theoretical indication that the melting of 
a uniform double-stranded polymer should be discontin- 
uous due to excluded volume effects [23[ . The discontinu- 
ity may be masked in experiments because of the inherent 
inhomogeneity of DNA [25|, [3(|, or by finite-size effects. 
The rigidity of DNA should play an important role in 



APPENDIX B: SLOPE OF THE PHASE 
BOUNDARY 

When the joint stiffness Jj vanishes, the partition 
function matrices and Tb in Eqs. (| 16[) and (JTTj) reduce 
to real-valued functions. As discussed in the Sec. IIII| 
two conditions have to be met at the phase transition, 
TsTi? = 1 and d z TB = oo, where in the former equation 
all multiplicative factors from the joints are absorbed in 
either of the two partition functions. Together, these two 
conditions set the value of /x = log(z) and en : = log(iu) 
along the phase boundary. All manipulations below are 
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then performed as the boundary point is changed by 
varying the force /. 

Noting that, the first condition is equivalent to 



iogr B (/x,/) + iogr fl ( M ,/,e R ) = o 



(Bl) 



its variations are obtained, by taking one total derivative 
with respect to /, as 



=8 f logT B + 8,, logT B d f n+ 

d f log T R + dfj, log T R dffi + d ER log T B d f e R 
=(L B ) + {N B )d f iJ+ 
(L R } + {N R )(d f (M + d f e R ). 

From the second condition we get 

dfdpTB _ (L B N B ) 



(B2) 



dadu^B 



(N B ) 



(B3) 



But because the condition d^Ts — oo of infinite bubble 
length is the same as log Tb = oo, one can equivalcntly 
express 



d f fx\ 



_dfdfJog£B_ 



(LbNb)c 

(N b )c 



(B4) 



where the subscript 'c' indicates a cumulant in place of 
a moment. From combining both expressions it follows 
that 



d f n 



d„r B =oo 



(Lb) 
(NbY 



(B5) 



Plugging this result into Eq. (|B2[) , one obtains 



d f e R 



(L, 



(Lr) 

(n b ) (N R y 



(B6) 



Note that at zero force all averages with only one L R or 
L R are zero by symmetry. Taking another total deriva- 
tive of Eq. (|B2[) . and dropping the terms that vanish for 
this reason, one finds that at / = 



0=(L|) c + (7V B )9 i 2 M + 

(L R ) c +(N R )(d 2 f n + d 2 f e R ). 



(B7) 



Taking another derivative of /i in Eq. (|B5|) one finds 
(for / = 0) 



d u r B =<x> 



(Nb)' 



(B8) 



Combining Eqs. (|B7|1 and (|B8[) the desired slope at / = 
is obtained as 



^2 _ ( L b)c _ ( L r)c 

f R (Nb) (N r ) ■ 



(B9) 
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